Phase separation in the 2D Hubbard model: 
a fixed-node quantum Monte Carlo study 
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Fixed-node Green's function Monte Carlo calculations have been performed for very large 16 x 16 
2D Hubbard lattices, large interaction strengths (7 = 10,20, and 40, and many (15 ~ 20) densities 
between empty and half filling. The nodes were fixed by a simple Slater- Gutz wilier trial wavefunc- 
tion. For each value of U we obtained a sequence of ground-state energies which is consistent with 
the possibility of a phase separation close to half-filling, with a hole density in the hole-rich phase 
which is a decreasing function of U. The energies suffer, however, from a fixed-node bias: more 
accurate nodes are needed to confirm this picture. Our extensive numerical results and their test 
against size, shell, shape and boundary condition effects also suggest that phase separation is quite 
a delicate issue, on which simulations based on smaller lattices than considered here are unlikely to 
give reliable predictions. 

PACS numbers: 71.10.Fd, 71.45.Lr, 74.20.-z 



Strongly correlated electrons and holes are expected 
to play a key role in the high-T c superconductors. Their 
possible instability towards phase separation (PS), ini- 
tially believed to inhibit superconductivity, is attracting 
a lot of interest since a few different authors have 
pointed out that such a tendency may in fact be inti- 
mately related to the high-T c superconductivity Long- 
range repulsive interactions may turn the PS instability 
into an incommensurate charge-density- wave (ICDW) in- 
stability, and the very existence of a quantum critical 
point associated to it may be a crucial ingredient of the 
superconducting transition Qj. PS and/or ICDW insta- 
bilities are related to a substantial reduction of the ki- 
netic energy, which otherwise tends to stabilize uniformly 
distributed states; such a reduction is typical of strongly 
correlated electrons, both in real and model systems. 

PS has been experimentally observed in h^CuOi-^s 
|^],[|, where the oxygen ions can move: in the doping in- 
terval 0.01 < S < 0.06 the compound separates into a 
nearly stoichiometric antiferromagnetic phase and a su- 
perconducting oxygen-rich phase. In generic compounds, 
where charged ions cannot move, the possibility of a 
macroscopic PS is spoiled by the long-range Coulomb 
repulsion, and should lead to an incommensurate CDW 
instability (?]]; here the identification of charge inhomo- 
geneities with spoiled PS is less straightforward ||. On 
the theoretical side, evidence for PS has been suggested 
for various models of strongly correlated electrons, as 
the t — J model (|], the three-band Hubbard model, the 
Hubbard-Holstein model and the Kondo model (see e.g. 
Ref. j| and references therein). 

Despite intensive studies, even for simple models there 
is no general agreement on the PS boundary: for the very 
popular t — J model, PS is fully established only at large 



J, but at small J (which unfortunately happens to be 
the physically relevant case) theoretical and numerical 
results are quite controversial. Emery et al.'s Q theory 
that PS occurs at any value of J in the t — J model is 
confirmed by a recent numerical study by Hellberg and 
Manousakis M], but is in contrast with Dagotto et al.'s 
H exact numerical results on small clusters, suggesting 
no tendency toward PS for both the Hubbard model and 
the t — J model below a critical value J < J c ~ t, and 
with Shih et al.'s JTT) numerical results. We also men- 
tion the recent suggestion by Gang Su |L2| , according to 
which the Hubbard model does not show PS at any value 
of U/t for any finite temperature, although it does not 
apply to ground-state properties. 

PS is a thcrmodinamic instability associated to the vi- 
olation, in a given density range n\ < n < U2, of the 
stability condition \ l — d 2 £/dn 2 > 0, which requires 
the energy density £ of an infinite electronic system to be 
a convex function of the electron density n. The system 
will therefore separate into two subsystems with electron 
densities n\ and ni- For the two-dimensional t-J and 
Hubbard models, PS, if any, is expected to occur in a 
density range close to half filling (n ~ 1), and to yield a 
hole-rich phase with density n\ < 1 and a hole-free phase 
with density rz.2 = 1 |». In a truly infinite system such 
a PS would be associated with a vanishing inverse com- 
pressiblity \~ l m t nc whole density range n\ < n < 112] 
in a finite system may even become negative, because 
of surface effects. So for finite systems it's preferable to 
pinpoint the PS using a Maxwell's construction (origi- 
nally suggested by Emery, Kivelson and Lin ||, see also 
below). But even such a procedure can give reliable re- 
sults only for medium-large finite systems; really small 
systems (for which most numerical results have been up 
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to now available) can attain so few and coarse densities, 
and suffer from so large finite-size errors, that their pre- 
dictions on the relevant trends remains largely inconclu- 
sive. 

Under these circumstances the availability of the fixed- 
node Green's function Monte Carlo (FNMC), a new and 
powerful numerical technique p3[ which allows the study 
of (previously unfeasible) large lattice- fermion systems, 
provides us with a powerful tool to further investigate 
the Hubbard model. Whether the 2D Hubbard hamilto- 
nian, a prototype for interacting electrons with no long- 
range repulsion, shows any instability towards PS, is a 
very interesting open question. A numerical study may 
also shed some indirect light on two related issues: the 
t — J model in the physical region of small J, and the 
adequacy of the one-band Hubbard hamiltonian to catch 
an essential aspect of high-T c superconductors. 

To evaluate the ground-state energy of the Hubbard 
hamiltonian 
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we thus implemented the FNMC method recently pro- 
posed for lattice fermions by van Bemmel et al. Jl3| . p~4[ . 
which has been used by Boninsegni for frustrated Heisen- 
berg systems |l5[] and by Gunnarsson et al. for orbitally- 
degenerate Hubbard models ]l6| ]. 

The Green's function Monte Carlo, after a sufficiently 
long imaginary time, projects out the ground-state com- 
ponent of any initial wavefunction; apart from transient 
estimates, which for large systems appear to be haz- 
ardous unless the initial variational wavefunction is suf- 
ficiently close to the exact one, this method is therefore 
not directly usable for fermions in our Hubbard model 
(as well as any other model whose Green's function is 
not positive everywhere). The FNMC |l3|]l4fl replaces 
the true hamiltonian by an effective hamiltonian which 
confines the Monte Carlo random walk within a single 
nodal region (a region of the configuration space where 
the guiding wavefunction never changes sign), and, in 
analogy with the continuum case |r^ , |l8|] , it provides an 
upper bound for the true ground-state energy ju) . 

We also implemented the tecnique proposed in Ref. 
|l9f , which allows us to reproduce with a relatively small 
fixed number (f 00 ~ 200) of walkers equally accurate re- 
sults as those obtained by means of standard MC runs of 
more than 2000 walkers. 

The variational wavefunctions we use to guide the ran- 
dom walks and to fix the nodes are the product of a 
Gutzwillcr factor and two Slater determinants of single- 
particle, mean-field wavefunctions for up- and down- 
spin electrons. The optimal Gutzwiller parameter and 
mean-field wavefunctions (whose only parameter is the 
staggered magnetization) were preliminarly obtained, for 
each U and density, by variational Monte Carlo (VMC) 
runs. 



A few representative variational and FNMC energies 
are shown in Table | for the 4x4 Hubbbard lattice, for 
which exact results [20] are available. As expected, the 
VMC energy is always above the FNMC energy, which 
for these coupling strengths is slightly (~ 3%) above the 
exact energy. For comparison we show the Constrained- 
Path Monte Carlo (CPMC) energies of Zhang et al. [§l), 
which also include a larger 16 x 16 lattice (last row). Es- 
pecially at large ?7's our results appear of comparable 
quality as theirs. As far as the 4x4 results are con- 
cerned, we notice that for N e — 10, which corresponds 
to a closed-shell configuration, both FNMC and CPMC 
are much closer to the exact energy than for N e = 14, 
which corresponds to an open-shell configuration. This 
could be a serious problem when numerically studying 
the behavior of the energy as a function of the density; 
the results presented here fortunately show that, for lat- 
tices larger than 12x12, the shell effects become almost 
irrelevant. 

To study the energy as a function of the electron den- 
sity we have first tried out the less usual way of vary- 
ing the density suggested by Ref. to avoid spurious 
Fermi-surface shape effects (keep the number of electrons 
N e fixed while the size of the underlying lattice is being 
varied), but discovered that either the number of elec- 
trons is really small (e.g. N e = 16), and then artificial 
changes in the convexity of the curve may occur, or the 
system is large enough (e.g. 12x12 lattices or larger), and 
then it doesn't matter how the density is being varied. 
So for our systematic study (many densities and three 
U values) we stick to a large 16 x 16 lattice (N s = 256 
sites), and vary the number of electrons N e to yield elec- 
tronic densities n — N e /N s ranging from empty n = to 
half filling n = 1. In the first panel of Fig. [l] we show 
the electronic ground-state energy per site, obtained by 
FNMC runs as a function of the density |22j] . Energies 
are in units of the hopping parameter t throughout this 
paper; the statistical errors are smaller than the marker 
size, and thus are not visible here. The calculated points 
are shown as full markers for closed shells, and as empty 
markers for open shells. From Fig. |] (see caption) it 
appears that the open-shell error, significant for a small 
4x4 lattice (see Table Q), becomes of the order of the sta- 
tistical error (and thus negligible) for our large lattices 
II- 

At all densities our three sets of data for U = 10 
(lower), 20 (middle), and 40 (upper curve) are brack- 
eted by the noninteracting unpolarized energy and the 
fully spin-polarized energy (both dashed in Fig. ^), and 
display a smooth and reasonable behavior. To evaluate 
the absolute accuracy of our results, we can rely on two 
exact limits: the low-density (n ~ 0) regime, where we 
expect £ — — 4n, and the half filled case (n = l),for which 
the strong-coupling expansion provides the correct large 
U behavior: to leading order in t/U, the model maps 
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onto an Heisenberg model, whose ground-state energy 
has been evaluated with great accuracy (l^j2^]. We can 
also consider the next correction term 34.6i 4 /£/ 3 j2j|. At 
low density our results are essentially exact; at half filling 
our error is small (~3%) for U = 10 but (as already seen 
in Table |) it tends to grow with U: ~ 9% for U = 20 
and ~11% for £7 = 40. We have made sure (see markers 
other than dots in Fig. |l| and footnote p2|) that such an 
energy discrepancy is not due to finite-size, shape, open- 
shell, and boundary condition effects; as far as systematic 
errors are concerned, we are thus left with the fixed-node 
approximation: as U grows, more flexible trial wavefunc- 
tions than adopted here are required to obtain accurate 
nodes [^6| . 

Keeping in mind the virtues and limitations of our nu- 
merical study, we can now turn to PS in the Hubbard 
model, ft has been shown in Ref. that the Maxwell 
construction is equivalent to study, as a function of the 
hole density x = l — n, the quantity e(x) = [eh{x)— eu]/x, 
i.e. the energy per hole eh(x) measured relative to its 
value at half filling en = eh(x = 0). For an infinite sys- 
tem, if the inverse compressibility vanishes between 
some critical density n c < 1 and half filling n = 1, then 
for < x < x c the function e(x) is a constant, and the 
fingerprint of a PS is thus a horizontal plot of e(x) be- 
low x c . For a finite system, instead, the PS fingerprint 
is a minimum of e{x) at x — x c P]. In some sense, e{x) 
works like a magnifying lens of PS. It should be stressed 
that in a consistent definition of e(x) the half- filling en- 
ergy en must be obtained as eh(x = 0) from the same 
calculation as any other eh{x 0) (in this work, from 
our FNMC). If that's not the case, then e(x) may tend 
to diverge near x = 0, with the danger of artificially cre- 
ating, rather than magnifying, the occurrence of PS. In 
the three right panels of Fig. [y we find plots of e{x) for 
U = 10, 20, and 40; these values, as well as the asso- 
ciated error bars, are obtained from those of the first 
panel (original FNMC energies and tiny error bars). De- 
spite the error bars, a common trend is evident for all the 
calculated coupling strenghts: e{x) has a positive slope 
for large hole densities, far from half- filling, but then it 
clearly changes slope below some small critical density 
x c . Such a minimum in e(x) implies that, at least for the 
FNMC effective hamiltonian determined by our choice of 
wavefunction, PS occurs below x — x c |p7j . Although a 
finer grid of hole densities would be required to locate 
with high precision the critical density function 
of U, we already see that x c decreases as U is increased; 
this qualitatively agrees with the original predictions || 
and with some previous calculations on the t — J model 
at corresponding values of J = At 2 /U JlCfl . 

In summary, our extensive FNMC numerical simula- 
tions of the Hubbard model for 16 x 16 two-dimensional 
lattices suggest PS for U t. If confirmed by further 
fixed-node simulations based on different nodes [ p6[ (and 
possibly even larger lattices pTJ), this result would imply 



that the t—J model is also likely to show PS in the physi- 
cally relevant regime J < 0.4, and that even a single-band 
Hubbard model is sufficient to reproduce this physical 
tendency of high-T c superconductors. 
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TABLE I. Ground-state energy per site (in units of the hopping parameter i) for a 4x4 Hubbard lattice and various values 
of U. N e is the number of electrons and n is the corresponding average density. VMC: variational Monte Carlo, this work; 
FNMC: Fixed-Node Green's function Monte Carlo, this work; CPMC: Constrained-Path Monte Carlo, Ref. [21]; EXACT: exact 
diagonalization results, Ref. [20]. (see text) 
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FIG. 1. The first panel, to the very left, shows the ground-state energy per site (in units of the hopping parameter t) as a 
function of the electronic density, for a 2D Hubbard lattice of N s — 16x16 = 256 sites with U — 10 (lower), 20 (middle), and 40 
(upper data). Errors are smaller than the marker size. Full markers correspond to closed shells and empty markers correspond 
to open shells. The dashed curves correspond to two (analytically given) U = results: the fully spin-polarized case (upper 
curve), whose total energy per site is symmetric with respect to quarter filling, and the unpolarized case (lower curve), whose 
total energy per site is symmetric with respect to half filling. Triangles (corresponding to a smaller 12 x 12 lattice and U = 40) 
and crosses (lly/2 x ll\/2 lattice and U = 10) are shown for comparison (see text). The second, third and fourth panel to the 
right contain plots of e(x) vs. x (see text) for U = 10, 20, and 40, respectively. The data markers have the same meaning as 
in the first panel; obviously at small x the error bar associated to e(x), Ae(x) = [Aej(i) + Aeh(x = 0)]/x, becomes significant 
even if the statistical FNMC error Aeh(i) is tiny (see text). 
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